/*------------------------------------------------------------------------
 *  Copyright (C) 2010  Luis M. de la Cruz
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

namespace Tuna {
/*!
 *  Write ScalarFields(1D, 2D, 3D) in ASCII format: x y z field for Parallel
 */
template<class T_sfield>
void InOut::writeToFileP(T_sfield& scf, int number, string filename,
			double dx, double Lx, 
			double dy = 1, double Ly = 0, 
			double dz = 1, double Lz = 0,
			int w = 5) 
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if( scf.rank() == 1 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << Lx + dx * i << "   " << scf(i) << endl;	
    }

    if ( scf.rank() == 2 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i) {	
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)	
	      file << Lx + dx * i << "   "		
		   << Ly + dy * j << "   " 
		   << scf(i,j) << endl;
	    file << endl;
	}
    }

    if (scf.rank() == 3 ) {
	for(int k = 0; k <= scf.ubound(thirdDim); ++k)
	    for(int i = 0; i <= scf.ubound(firstDim); ++i)
		for(int j = 0; j <= scf.ubound(secondDim); ++j)
		    file << Lx + dx * i << "   " 
			 << Ly + dy * j << "   "
			 << Lz + dz * k << "   " 
			 << scf(i,j,k) << endl;
    }

    file.close();
}


// Write ScalarFields(1D, 2D and 3D) in OpenDX format (Parallel)
///
/// This function must write scalar fields in 1 2 and 3 dimensions using
/// the data model of OpenDX for regular and rectangular meshes.
///
template<class T_sfield>
void InOut::writeToFile_DXP(T_sfield& scf, int number, string filename,
			    double dx = 1, double Lx = 0,
			    double dy = 1, double Ly = 0,
			    double dz = 1, double Lz = 0,
			    int w = 5) 
{   
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if( scf.rank() == 1 ) {
	file << "object 1 class array type float rank 1 shape 2 items "
	     << scf.size() << " data follows\n";
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << dx * i << "   " << scf(i) << endl;	
	file << "object 2 class array type int rank 1 shape 2 items "
	     << scf.size() - 1 << " data follows\n";
	for(int i = 0; i < scf.ubound(firstDim); ++i)
	    file << i << "   " << i + 1 << endl;
	file << "attribute \"element type\" string \"lines\"\n"
	     << "attribute \"ref\" string \"positions\"\n";
	file << "object 3 class array type float rank 0 items "
	     << scf.size() << " data follows\n";
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << scf(i) << "\t";
	file << endl;

	file << "attribute \"dep\" string \"positions\"\n";
	file << "object \"irregular positions irregular connections\" class field\n"
	     << "component \"positions\" value 1\n" 
	     << "component \"connections\" value 2\n" 
	     << "component \"data\" value 3\n" 
	     << "end\n" ;
    }
    if ( scf.rank() == 2 ) {
	int Nx = scf.extent(firstDim);
	int Ny = scf.extent(secondDim);
	int count;
	file << "object 1 class gridpositions counts " 
	     << Nx << " " << Ny << endl;
	file << "origin \t " << Lx << " \t " << Ly << " \n";
	file << "delta \t " << dx << " \t 0 \n";
	file << "delta \t 0 \t " << dy << "\n";
	file << "object 2 class gridconnections counts " 
	     << Nx << " " << Ny << endl;
	file << "attribute \"element type\" string \"quads\" \n";
	file << "attribute \"ref\" string \"positions\" \n";
	file << "object 3 class array type float rank 0 items " 
	     << Nx * Ny << " data follows\n";
	file.precision(10);
	file.setf(ios_base::fixed,ios_base::floatfield);	
	for(int i = 0; i <= scf.ubound(firstDim); ++i) 
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)
		file << scf(i,j) << "\t";
	file << "\nattribute \"dep\" string \"positions\" \n";
	file<<"object \"regular positions regular connections\" class field\n";
	file << "component \"positions\" value 1 \n";
	file << "component \"connections\" value 2 \n";
	file << "component \"data\" value 3 \n";
	file << "end\n";
    }
    if (scf.rank() == 3 ) {
	int Nx = scf.extent(firstDim);
	int Ny = scf.extent(secondDim);
	int Nz = scf.extent(thirdDim);
	int count;
	file << "object 1 class gridpositions counts " 
	     << Nx << " " << Ny << " " << Nz << endl;
	file << "origin \t 0 \t 0 \t 0 \n";
	file << "origin \t " << Lx << " \t " << Ly << " \t " << Lz << " \n";
	file << "delta \t " << dx << " \t 0 \t 0 \n";
	file << "delta \t 0 \t " << dy << " \t 0 \n";
	file << "delta \t 0 \t 0 \t " << dz << "\n";
	file << "object 2 class gridconnections counts " 
	     << Nx << " " << Ny << " " << Nz << endl;
	file << "attribute \"element type\" string \"cubes\" \n";
	file << "attribute \"ref\" string \"positions\" \n";
	file << "object 3 class array type float rank 0 items " 
	     << Nx * Ny * Nz << " data follows\n";
	file.precision(10);
	file.setf(ios_base::fixed,ios_base::floatfield);
	
	for(int i = 0; i <= scf.ubound(firstDim); ++i) 
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)
		for(int k = 0; k <= scf.ubound(thirdDim); ++k) 
		    file << scf(i,j,k) << "\t";
	file << "\nattribute \"dep\" string \"positions\" \n";
	file<<"object \"regular positions regular connections\" class field\n";
	file << "component \"positions\" value 1 \n";
	file << "component \"connections\" value 2 \n";
	file << "component \"data\" value 3 \n";
	file << "end\n";
    }
    file.close();
}

  
/*!
 *  Write ScalarFields(1D, 2D, 3D) in ASCII format: x y z field
 */
template<class T_sfield>
void InOut::writeToFile(T_sfield& scf, int number, string filename,
			 double dx, double dy = 1, double dz = 1,
			 int w = 5) 
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if( scf.rank() == 1 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << dx * i << "   " << scf(i) << endl;	
    }

    if ( scf.rank() == 2 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i) {	
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)	
		file << dx * i << "   "		
		     << dy * j << "   " 
		     << scf(i,j) << endl;
		file << endl;
		}
    }

    if (scf.rank() == 3 ) {
	for(int k = 0; k <= scf.ubound(thirdDim); ++k)
	    for(int i = 0; i <= scf.ubound(firstDim); ++i)
		for(int j = 0; j <= scf.ubound(secondDim); ++j)
		    file << dx * i << "   " 
			 << dy * j << "   "
			 << dz * k << "   " 
			 << scf(i,j,k) << endl;
    }

    file.close();
}

// Read files written by the above writeFile() function.
template<class T_sfield>
void InOut::readFromFile(T_sfield& scf, int number, string filename, 
			 double dx, double dy = 1, double dz = 1,
			 int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

// ASCII InOut
    ifstream file (filename.c_str());

// BINARY InOut
//    ofstream file(filename.c_str(), ios::out | ios::app | ios::binary);
	
    if (!file) {
	cout << "\n +---------------+\n"
	     << " * TUNA-CFD Error: !Can't open file \"" 
	     << filename << "\" * \n +---------------+\n\n";
    }
    
    if( scf.rank() == 1 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file >> dx >> scf(i) ;	
    }
    

    if ( scf.rank() == 2 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i) {
	    for(int j = 0; j <= scf.ubound(secondDim); ++j) 
		file >> dx >> dy >> scf(i,j);
	}
    }

    if (scf.rank() == 3 ) {
	for(int k = 0; k <= scf.ubound(thirdDim); ++k)
	    for(int i = 0; i <= scf.ubound(firstDim); ++i) 
		for(int j = 0; j <= scf.ubound(secondDim); ++j) {
		    file >> dx >> dy >> dz >> scf(i,j,k);
		}
    }

    file.close();
}
    
// Write ScalarFields(1D, 2D, 3D) in ASCII format: field
// This function doesn't store dx, dy and dz. 
template<class T_sfield>
void InOut::writeToFile(T_sfield& scf, int number, string filename, int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;


// ASCII InOut
    ofstream file (filename.c_str());

// BINARY InOut
//    ofstream file(filename.c_str(), ios::out | ios::app | ios::binary);
	
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if( scf.rank() == 1 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << scf(i) << "   ";	
    }
    

    if ( scf.rank() == 2 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i) {
	    for(int j = 0; j <= scf.ubound(secondDim); ++j) 
		file << scf(i,j) << "   ";
	}
    }

    if (scf.rank() == 3 ) {
	for(int k = 0; k <= scf.ubound(thirdDim); ++k)
	    for(int i = 0; i <= scf.ubound(firstDim); ++i) 
		for(int j = 0; j <= scf.ubound(secondDim); ++j) {
		    file << scf(i,j,k) << "   ";
		}
    }
    file.close();
}


// Read files written by the above writeFile() function.
// This function doesn't read dx, dy and dz.
template<class T_sfield>
void InOut::readFromFile(T_sfield& scf, int number, string filename, int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

// ASCII InOut
    ifstream file (filename.c_str());

// BINARY InOut
//    ofstream file(filename.c_str(), ios::out | ios::app | ios::binary);
	
    if (!file) {
	cout << "\n +---------------+\n"
	     << " * TUNA-CFD Error: !Can't open file \"" 
	     << filename << "\" * \n +---------------+\n\n";
    }
    
    if( scf.rank() == 1 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file >> scf(i) ;	
    }
    

    if ( scf.rank() == 2 ) {
	for(int i = 0; i <= scf.ubound(firstDim); ++i) {
	    for(int j = 0; j <= scf.ubound(secondDim); ++j) 
		file >> scf(i,j);
	}
    }

    if (scf.rank() == 3 ) {
	for(int k = 0; k <= scf.ubound(thirdDim); ++k)
	    for(int i = 0; i <= scf.ubound(firstDim); ++i) 
		for(int j = 0; j <= scf.ubound(secondDim); ++j) {
		    file >> scf(i,j,k);
		}
    }
    file.close();
}
    

// Write ScalarFields(1D, 2D and 3D) in OpenDX format
/// 
/// This function must write scalar fields in 1 2 and 3 dimensions using
/// the data model of OpenDX for regular and rectangular meshes.
///
template<class T_sfield>
void InOut::writeToFile_DX(T_sfield& scf, int number, string filename,
			    double dx = 1, double dy = 1, double dz = 1,
			    int w = 5) 
{   
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if( scf.rank() == 1 ) {
	file << "object 1 class array type float rank 1 shape 2 items "
	     << scf.size() << " data follows\n";
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << dx * i << "   " << scf(i) << endl;	
	file << "object 2 class array type int rank 1 shape 2 items "
	     << scf.size() - 1 << " data follows\n";
	for(int i = 0; i < scf.ubound(firstDim); ++i)
	    file << i << "   " << i + 1 << endl;
	file << "attribute \"element type\" string \"lines\"\n"
	     << "attribute \"ref\" string \"positions\"\n";
	file << "object 3 class array type float rank 0 items "
	     << scf.size() << " data follows\n";
	for(int i = 0; i <= scf.ubound(firstDim); ++i)
	    file << scf(i) << "\t";
	file << endl;

	file << "attribute \"dep\" string \"positions\"\n";
	file << "object \"irregular positions irregular connections\" class field\n"
	     << "component \"positions\" value 1\n" 
	     << "component \"connections\" value 2\n" 
	     << "component \"data\" value 3\n" 
	     << "end\n" ;
    }
    if ( scf.rank() == 2 ) {
	int Nx = scf.extent(firstDim);
	int Ny = scf.extent(secondDim);
	int count;
	file << "object 1 class gridpositions counts " 
	     << Nx << " " << Ny << endl;
	file << "origin \t 0 \t 0 \n";
	file << "delta \t " << dx << " \t 0 \n";
	file << "delta \t 0 \t " << dy << "\n";
	file << "object 2 class gridconnections counts " 
	     << Nx << " " << Ny << endl;
	file << "attribute \"element type\" string \"quads\" \n";
	file << "attribute \"ref\" string \"positions\" \n";
	file << "object 3 class array type float rank 0 items " 
	     << Nx * Ny << " data follows\n";
	file.precision(10);
	file.setf(ios_base::fixed,ios_base::floatfield);	
	for(int i = 0; i <= scf.ubound(firstDim); ++i) 
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)
		file << scf(i,j) << "\t";
	file << "\nattribute \"dep\" string \"positions\" \n";
	file<<"object \"regular positions regular connections\" class field\n";
	file << "component \"positions\" value 1 \n";
	file << "component \"connections\" value 2 \n";
	file << "component \"data\" value 3 \n";
	file << "end\n";
    }
    if (scf.rank() == 3 ) {
	int Nx = scf.extent(firstDim);
	int Ny = scf.extent(secondDim);
	int Nz = scf.extent(thirdDim);
	int count;
	file << "object 1 class gridpositions counts " 
	     << Nx << " " << Ny << " " << Nz << endl;
	file << "origin \t 0 \t 0 \t 0 \n";
	file << "delta \t " << dx << " \t 0 \t 0 \n";
	file << "delta \t 0 \t " << dy << " \t 0 \n";
	file << "delta \t 0 \t 0 \t " << dz << "\n";
	file << "object 2 class gridconnections counts " 
	     << Nx << " " << Ny << " " << Nz << endl;
	file << "attribute \"element type\" string \"cubes\" \n";
	file << "attribute \"ref\" string \"positions\" \n";
	file << "object 3 class array type float rank 0 items " 
	     << Nx * Ny * Nz << " data follows\n";
	file.precision(10);
	file.setf(ios_base::fixed,ios_base::floatfield);
	
	for(int i = 0; i <= scf.ubound(firstDim); ++i) 
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)
		for(int k = 0; k <= scf.ubound(thirdDim); ++k) 
		    file << scf(i,j,k) << "\t";
	file << "\nattribute \"dep\" string \"positions\" \n";
	file<<"object \"regular positions regular connections\" class field\n";
	file << "component \"positions\" value 1 \n";
	file << "component \"connections\" value 2 \n";
	file << "component \"data\" value 3 \n";
	file << "end\n";
    }
    file.close();
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Write ScalarFields(2D) that are components of 2D velocity in OpenDX format (Parallel)
template<class T_sfield>
void InOut::writeToFile_DXP(T_sfield& scf1, T_sfield& scf2,
			   int number, string filename, 
			   double dx = 1, double Lx = 0,
			   double dy = 1, double Ly = 0,
			   int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');   
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }

    int Nx = scf1.extent(firstDim);
    int Ny = scf1.extent(secondDim);
    int count;
    file << "object 1 class gridpositions counts " 
	 << Nx << " " << Ny << endl;
    file << "origin \t " << Lx << " \t " << Ly << " \n";
    file << "delta \t " << dx << " \t 0 \n";
    file << "delta \t 0 \t " << dy << " \n";
    file << "object 2 class gridconnections counts " 
	 << Nx << " " << Ny  << endl;
    file << "attribute \"element type\" string \"quads\" \n";
    file << "attribute \"ref\" string \"positions\" \n";
    file << "object 3 class array type float rank 1 shape 2 items " 
	 << Nx * Ny << " data follows\n";

    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);    
    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
      for(int j = 0; j <= scf1.ubound(secondDim); ++j) {
	file << scf1(i,j) << "\t" 
	     << scf2(i,j) << "\t";
	file << endl;
      }

    file << "\nattribute \"dep\" string \"positions\" \n";
    file << "object \"regular positions regular connections\" class field\n";
    file << "component \"positions\" value 1 \n";
    file << "component \"connections\" value 2 \n";
    file << "component \"data\" value 3 \n";
    file << "end\n";  

    file.close();
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

// Write ScalarFields(2D) that are components of 2D velocity in OpenDX format 
template<class T_sfield>
void InOut::writeToFile_DX(T_sfield& scf1, T_sfield& scf2,
			   int number, string filename, 
			   double dx = 1, double dy = 1, 
			   int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');   
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }

    int Nx = scf1.extent(firstDim);
    int Ny = scf1.extent(secondDim);
    int count;
    file << "object 1 class gridpositions counts " 
	 << Nx << " " << Ny << endl;
    file << "origin \t 0 \t 0 \n";
    file << "delta \t " << dx << " \t 0 \n";
    file << "delta \t 0 \t " << dy << " \n";
    file << "object 2 class gridconnections counts " 
	 << Nx << " " << Ny  << endl;
    file << "attribute \"element type\" string \"quads\" \n";
    file << "attribute \"ref\" string \"positions\" \n";
    file << "object 3 class array type float rank 1 shape 2 items " 
	 << Nx * Ny << " data follows\n";
    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);
    
    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j) {
	    file << scf1(i,j) << "\t" 
		 << scf2(i,j) << "\t";
	    file << endl;
	}

    file << "\nattribute \"dep\" string \"positions\" \n";
    file << "object \"regular positions regular connections\" class field\n";
    file << "component \"positions\" value 1 \n";
    file << "component \"connections\" value 2 \n";
    file << "component \"data\" value 3 \n";
    file << "end\n";  

    file.close();
}


// Write ScalarFields(3D) that are components of 3D velocity in OpenDX format 
template<class T_sfield>
void InOut::writeToFile_DX(T_sfield& scf1, T_sfield& scf2, T_sfield& scf3,
			   int number, string filename, 
			   double dx = 1, double dy = 1, double dz = 1,
			   int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');   
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }

    int Nx = scf1.extent(firstDim);
    int Ny = scf1.extent(secondDim);
    int Nz = scf1.extent(thirdDim);
    int count;
    file << "object 1 class gridpositions counts " 
	 << Nx << " " << Ny << " " << Nz << endl;
    file << "origin \t 0 \t 0 \t 0 \n";
    file << "delta \t " << dx << " \t 0 \t 0 \n";
    file << "delta \t 0 \t " << dy << " \t 0 \n";
    file << "delta \t 0 \t 0 \t " << dz << "\n";
    file << "object 2 class gridconnections counts " 
	 << Nx << " " << Ny << " " << Nz << endl;
    file << "attribute \"element type\" string \"cubes\" \n";
    file << "attribute \"ref\" string \"positions\" \n";
    file << "object 3 class array type float rank 1 shape 3 items " 
	 << Nx * Ny * Nz << " data follows\n";
    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);
    
    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j) {
	    for(int k = 0; k <= scf1.ubound(thirdDim); ++k)
		file << scf1(i,j,k) << "\t" 
		     << scf2(i,j,k) << "\t"	
		     << scf3(i,j,k) << "\t" ;	
	    file << endl;
	}

    file << "\nattribute \"dep\" string \"positions\" \n";
    file << "object \"regular positions regular connections\" class field\n";
    file << "component \"positions\" value 1 \n";
    file << "component \"connections\" value 2 \n";
    file << "component \"data\" value 3 \n";
    file << "end\n";  

    file.close();
}


// Write ScalarField(1D, 2D y 3D) in OpenDX format for NonUniform Meshes 
template<class T_sfield, class T_mesh>
void InOut::writeToFile_DX(T_sfield& scf1, 
			   int number, string filename, T_mesh& m, int w = 5) {
    vector<double> x = m.getX(), y = m.getY(), z = m.getZ();
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
  
    if ( scf1.rank() == 1 ) {
      for(int i = 0; i <= scf1.ubound(firstDim); ++i)
	file << x[i] << "\t" << scf1(i) << endl;
	
    }  
    if ( scf1.rank() == 2 ) {
      for(int j = 0; j <= scf1.ubound(secondDim); ++j) 
	for(int i = 0; i <= scf1.ubound(firstDim); ++i)
	  file << x[i] << "\t" << y[j] << "\t" << scf1(i,j) << endl; 
    }
    if ( scf1.rank() == 3 ) {
      for(int k = 0; k <= scf1.ubound(thirdDim); ++k) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j) 
	  for(int i = 0; i <= scf1.ubound(firstDim); ++i)
	    file << x[i] << "\t" << y[j] << "\t" << z[k] << "\t" << scf1(i,j,k) << endl; 
    }

    file.close();
}

// Write ScalarFields(2D) that are components of a 2D vector field in OpenDX 
// format for NonUniform Meshes
template<class T_sfield, class T_mesh>
void InOut::writeToFile_DX(T_sfield &scf1, T_sfield &scf2, 
			   int number, string filename, T_mesh &m, int w = 5) {
    vector<double> x = m.getX(), y = m.getY(), z = m.getZ();
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');    
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if ( scf1.rank() == 2 ) {
	for(int j = 0; j <= scf1.ubound(secondDim); ++j) 
	    for(int i = 0; i <= scf1.ubound(firstDim); ++i)
		file << x[i] << "\t" << y[j] << "\t" 
		     << scf1(i,j) << "\t"
		     << scf2(i,j) << endl;	
    }

    file.close();
}

// Write ScalarFields(3D) that are components of a 3D vector field in OpenDX 
// format for NonUniform Meshes
template<class T_sfield, class T_mesh>
void InOut::writeToFile_DX(T_sfield &scf1, T_sfield &scf2, T_sfield &scf3, 
			   int number, string filename, T_mesh &m, int w = 5) {
    vector<double> x = m.getX(), y = m.getY(), z = m.getZ();
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');    
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    if ( scf1.rank() == 3 ) {
      for(int k = 0; k <= scf1.ubound(thirdDim); ++k) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j) 
	  for(int i = 0; i <= scf1.ubound(firstDim); ++i)
	    file << x[i] << "\t" << y[j] << "\t" << z[k] << "\t" 
		 << scf1(i,j,k) << "\t"
		 << scf2(i,j,k) << "\t"
		 << scf3(i,j,k) << "\t"
		 << endl; 
    }

    file.close();
}

// Amiramesh file format
// This function writes amira mesh format for the 3D case
template<class T_sfield>
void InOut::writeToFile_AM(T_sfield& scf, int number, string filename,
			   double dx = 1, double dy = 1, double dz = 1,
			   int w = 5) 
{   
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ofstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }

    if (scf.rank() == 3 ) {
	int Nx = scf.extent(firstDim);
	int Ny = scf.extent(secondDim);
	int Nz = scf.extent(thirdDim);
	int count;
	file << "# AmiraMesh 3D ASCII 2.0 \n"
	     << "# CreationDate: Sun Set 30 2007 \n"
	     << "define Lattice " << Nx << " " << Ny << " " << Nz << endl
	     << "Parameters { \n"
	     << "    CoordType \"uniform\", \n"
	     << "    BoundingBox 0 1 0 1 0 1 \n"
	     << "}\n\n"
	     << "Lattice { float ScalarField } = @1 \n"
	     << "@1 \n";
	file.precision(10);
	file.setf(ios_base::fixed,ios_base::floatfield);
	for(int i = 0; i <= scf.ubound(firstDim); ++i) 
	    for(int j = 0; j <= scf.ubound(secondDim); ++j)
		for(int k = 0; k <= scf.ubound(thirdDim); ++k) 
		    file << scf(i,j,k) << "\t";
	file << "\n";
    }

    file.close();
}



template<class T_sfield>
void InOut::readFromFile_1(T_sfield& scf1, T_sfield& scf2, 
			    int number, string filename, int w = 5)
{
    ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;

    ifstream file (filename.c_str());
    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
///
/// Here, scf1 and scf2 are supposed to be of type ScalarField2D
///
    for(int j = 0; j <= scf1.ubound(secondDim); ++j) 
	for(int i = 0; i <= scf1.ubound(firstDim); ++i)
	    file >> scf1(i,j) >> scf2(i,j) ;
}

///
/// Funtion to read from OpenDX native files. This function just read
/// files written using printToFile_DXNative, and not any kind of 
/// DX files.
/// Last modification: Wed Sep  6 17:56:51 CDT 2006
///
///
template<class T_sfield>
void InOut::readFromFile_2(T_sfield& scf1, T_sfield& scf2, T_sfield& scf3, 
			   int number, string filename, int w = 5) 
{   
  ostringstream extnumber;
    extnumber.width(w);
    extnumber.fill('0');
    extnumber << number;
    string snumber = extnumber.str();
    filename += snumber;
    ifstream file (filename.c_str());

    if (!file) {
	cout << "\n !Can't open file \" " << filename << "\"" << endl;
    }
    
    char line[256];
/// 
/// by the time, just jump the first 7 lines
///	
    for(int ii = 0; ii < 9; ++ii)
	file.getline(line,256);
    
    file.precision(10);
    file.setf(ios_base::fixed,ios_base::floatfield);
///
/// Here, scf1, scf2 and scf3 are supposed to be of type ScalarField3D
///
    for(int i = 0; i <= scf1.ubound(firstDim); ++i) 
	for(int j = 0; j <= scf1.ubound(secondDim); ++j)
	    for(int k = 0; k <= scf1.ubound(thirdDim); ++k) {
		file >> scf1(i,j,k) >> scf2(i,j,k) >> scf3(i,j,k);
	    }
    file.close();
}




  
// Write out the Diagonal Matrices to a file. (1D)
ostream& operator<<(ostream& s, const DiagonalMatrix<float,1>& M)
{
    s << "\n a = " << M.a << "\n b = " << M.b << "\n c = " << M.c << endl;
    return s;
}
  
// Write out the Diagonal Matrices to a file. (2D)
ostream& operator<<(ostream& s, const DiagonalMatrix<float,2>& M)
{
    s << "\n a = " << M.a << "\n b = " << M.b << "\n c = " << M.c 
      << "\n d = " << M.d << "\n e = " << M.e << endl;    
    return s;
}

// Write out the Diagonal Matrices to a file. (3D)  
ostream& operator<<(ostream& s, const DiagonalMatrix<float,3>& M)
{
    s << "\n a = " << M.a << "\n b = " << M.b << "\n c = " << M.c  
      << "\n d = " << M.d << "\n e = " << M.e << "\n f = " << M.f
      << "\n g = " << M.g << endl;
    return s;
}


} // Tuna namespace
